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ABSTRACT 

We present an analysis of star formation and feedback recipes appropriate for galac- 
tic smoothed particle hydrodynamics simulations. Using an isolated Milky Way-like 
galaxy, we constrain these recipes based on well-established observational results. Our 
star formation recipe is based on that of Katz (1992) with the additional inclusion of 
physically motivated supernova feedback recipes. We propose a new feedback recipe 
in which type II supernovae are modelled using an analytical treatment of blastwaves. 
With this feedback mechanism and a tuning of other star formation parameters, the 
star formation in our isolated Milky Way-like galaxy follows the slope and normali- 
sation of the observed Schmidt law. In addition, we reproduce the low density cutoff 
and filamentary structure of star formation observed in disk galaxies. Our final recipe 
will enable better comparison of N-body simulations with observations. 

Key words: star formation - supernova feedback - smoothed particle hydrodynam- 
ics. 



1 INTRODUCTION 

Meaningful comparisons between observations and simula- 
tions of galaxies require that simulations include gas and 
stars. While dark matter controls the global dynamics of 
galaxies, it cannot be directly observed. Gas and stars pro- 
vide the photons that comprise observations. Unfortunately, 
for the foreseeable future, galaxy simulations will not be 
able to resolve individual stars, their individual explosions 
as supernovae, or the fine structure of gas clouds in the 
context of a galaxy. It is only possible to use a heuristic 
recipe to model the subresolution physics of star forma- 
tion and the feedback effect of supernova explosions. Sim- 
ulations that use the simplest methods of injecting super- 
nova energy into the gas surrounding stars have proven 
ineffective at producing realistic star formation feedback 
!Katil992lh The introduction of more clever sub resolution 
schemes for the distribut ion of supernova energy as a feed- 
back on star formation jThacker fc Couchmanl i200 and 
ISpringel fc Hernauistl <l2003fO has resulted in more realistic 
disk galaxies. 

The two methods currently employed to model gas 
physics in simulations are Eulerian grid codes and La- 
grangian particle codes. Eulerian codes (e.g. lCen fc OstrikeJ 
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Jl995l) . lKravtsolJl99^ . and lO'Shea et al.H2004M track the 
movement of gas around a fixed grid of cells. The method 
that we describe in this paper involves a Lagrangian treat- 
ment o f the gas called smoothed particle hydrodynamics 
(SPH) iMonagharj ri992h Resolution is achieved naturally 
as particles concentrate in dense regions of interest. To solve 
the equations of hydrodynamics, physical quantities are de- 
termined using spline kernel interpolation of neighbouring 
particles. 

Various methods have been employed to convert 
these smoothed gas particles into stars. Evidence of the 
improvemen t in the resolution of simulations is that 
early work JCen fc Ostrikerlll995t chan ged dense gas into 
"galaxy particles". More re c ent codes tygBgj_et_al J ^)9j1; 
iHultman fc Pharasvnl Il999t ISpringel fc Hernauistl [20031 
choose to package stars along with hot and cold gas into 
multiphase particles correspondi ng to the cold/warm an d 
hot phases of the ISM reported in iMcKee fc Ostrikerl Jl977t) . 
Physical properties like temperature and density of the mul- 
tiphase particles are assumed to be the average of the differ- 
ent components of ea ch particle. In one vers i on of m ultiphase 
particles ('explicit'), ISpringel fc Hernauistl (12003') allow gas 
particles to host stellar mass, spawning stars once the stel- 
lar component of the multiphase particle exceeds a minimum 
star particle mass. 

Using the SPH code GASOLINE JWadslev et alJl2004Tl . 
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w e choose to follow a different method that is outlined 
in IKatzl (Il992f) , one tha t uses a stochastic formu lation for 
star formation (see also Kawata fc Gibson! ll2003r). th e 'or- 
dinary' star formation in lSpringel fc Hernq uist (2003), and 
lOkamoto et al.l (|200^)). This allows stars to form immedi- 
ately and gas particles to maintain their own, unique charac- 
ter and temperature, hot, cold or warm. As supernovae are 
the direct result of star formation, we use them as the feed- 
back that limits star formation. The coarse mass and spatial 
resolution of current simulations limits models to heuristic 
descriptions of this phenomenon. The insufficient resolution 
and lack of multiple gas phases also means that the 10 51 ergs 
of energy that a supernova generates, if deposited as thermal 
energy, would be dissipated through radiative cooling pro- 
cesses before that thermal energy has any effect on the gas 
surrounding the s upernova explosion (e.g. iKatz et alJll99(i 
iBrook et al1l2004ft . 

Two metho ds have been emp l oyed t o use that energy 
in other ways. iNavarro fc Whita (I1993T) describe a kine- 
matic feedback mechanism that uses the supernova energy 
to provide an outwards velocity kick to all of the gas par- 
ticles surrounding a star particle in which a supernova (or 
more commonly a lar g e grou p of supernova) has exploded. 
ISpringel fc Hernauistl i2003fl use this idea to individually 
kick their multiphase particles in some direction, either ran- 
domly or perpendicular to their angular momentum vector, 
after a supernova has exploded. Even in such kinematic ap- 
proaches the feedback does not have a very strong effect, 
i.e. it is not efficient in driving winds from small galaxies, 
unles s the hydrodynamic forces a re temporarily turned off 
as in lSpringel fc Hernauistl <l2003l) . 

An alternative to these kinematic examples is to 
more effectively mimic the transfer of the kinetic en- 
ergy of supernova Shockwaves to thermal energ y in 
the interstellar medi u m (I SM). lYepes et al.l il997|) and 
iHultman fc PharasvrJ <ll999F) both utilise multiphase com- 
ponents and convert some multiple of the mass in stars 
undergoing type SNII explosions from a co ld gas phase 
to a separate hot gas phase. Alternatively, IPearce et alJ 
(1999) assigned gas particles to one of two phases in 
SPH simulations. The phase separation can limit the 
loss of thermal energy de posited from supernova feed- 
back jMarri fc White! l2003ft . In contrast, iGerritsenl dl997f) 
does not explicity use multiphase particles but turns off 
the radiative cooling of the gas particles immediately sur- 
roundin g a star particle in w h ich S N II have r e cently ex- 
ploded. IThacker fc Couchmanl fcOOOl) ; iBottemal (120031) ex- 
amined how the TGeTriteenT " ^9^) scheme works in a n 
isolated Milky Way galaxy. IThacker fc Couchmanl l|200llh 
ISommer-Larsen et al.l <[2003l)nGovernato et alJ j2006Tt used 
this recipe in cosmolog ical simulations of forming galaxies. 
iPelupessv et al.| (120041) explored how this reci pe worked in 
the en vironment of isolated disk dwarf galaxies. IBrook et alJ 
i2004Tl showed that the adiabatic feedback method produces 
more realistic feedback in simulations of isolated collapsing 
halos. 

We extend the e x plorat ion of the recipe started in 
IThacker fc Couchmanl |2000) for the case of an isolated 
Milky Way galaxy here with the goal of further refining the 
star formation and feedback algorithm. This paper specif- 
ic ally focuses on the star formation recipe first described 
in IKatzl il992l) paired with feedback that for a time cools 



only adiabat ically to prevent immediate r adiative losses as 
described in IThacker fc Couchmar] 12000^1 . Because of this 
focus, we are able to optimise the recipe, present a full pa- 
rameter study and include investigations of specific issues 
like the resolution dependence of star formation and feed- 
back. 

^describes the details of the recipe that we test to con- 
vert gas into stars and describes our feedback schemes. 
presents the tunable parameters that affect star forma- 
tion in our isolated model Milky Way (IMMW). presents 
the results of varying the criteria and parameters and deter- 
mines which choices work best at reproducing observations 
for the IMMW. [JS] discusses the relevance and plausibility 
of our recipe and possible future improvements. We present 
our final star formation and feedback recipe and conclude in 

m 



2 STAR FORMATION 

O ur star form ation algorithm is similar to the one proposed 
in IKatzl Jl992h and extended in IKatz et all Jl996l . hereafter 
KWH). First, we apply criteria to determine which gas par- 
ticles are eligible to form stars. We then determine which 
gas particles actually form stars probabilistically such that 
on average we reprodu ce a star format ion rate formula sim- 
ilar to a Schmidt law (ISchmiddll959T) . Those gas particles 
that actually form stars spawn a new star particle of a pre- 
determined mass, reducing their own mass accordingly. The 
new star particle is created with the same velocity, position, 
and metallicity as its parent gas particle. Star particles can 
add energy, mass and metals back to gas particles through 
feedback processes including type II and la supernova and 
stellar winds. 



2.1 Criteria 

The recipe proposed in IKatzl Jl992t) and KWH starts with 
an examination of every SPH gas particle in the simulation. 
The gas particle must satisfy 4 criteria before it is eligible 
for star formation: 



Is the particle denser than n m i n = 0.1 cm 
Is the particle in an overdense region? 
Is the particle part of a converging flow? 
Is the particle Jeans unstable (— > 



The density and overdensity criteria simply check that par- 
ticles fall above the limits set for the simulation. We choose 
the overdensity limit to be 55 p/p, which limits star forma- 
tion to virialised regions at early times in the Universe when 
the physical density everywhere is high and plays no role in 
simulations of isolated galaxies. It is a simple matter to ad- 
just these limits to match observed galaxy properties. In this 
paper we use n to represent the number density (cm -3 ) and 
p to represent the mass density (g cm -3 ). For gas particles 
npnriH = p, where p is the mean molecular weight and uih 
is t he mass of a Hydrogen atom. 

IKat j il992l) made the reasonable assumption that the 
gas forming a star should be in a collapsing region and so 
required that the gas particles be part of a converging flow. 
In our implementation of SPH every particle is assigned a 
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smoothing length, h, such that there are a fixed number 
of particles (neighbours), iV smoot h, within twice that length. 
We usually choose iV 8mo oth to be 32. Physical quantities are 
estimated using spline kernel interpolation. For example, the 
mass density, p, for particle i is given by 

JV 

P' 1 = ^ mjWQrj - Tj\, hj, hj) (1) 
j=i 

where m is the particle mass, N is the number of gas parti- 
cles, and W is the smoothing kernel, which we choose to have 
compact support, i.e. it goes to zero beyond 2h so the sum is 
really only over JV smoo th iMonaghanll992l : lHernauist fc Kat3 
1989). The divergence of the velocity field, V ■ v at the po- 
sition of gas particle i is given by 



V-V: 



1 2 m i( v i- v 0-ViW(|r 4 -rj|,/u,/»,-) (2) 



Pi 



j=i 



where v is the velocity. When V ■ v is negative, the criterion 
is satisfied, it is assumed that the gas particle is part of a 
collapsing flow, and the gas particle can form stars. 

The Jeans Criterion is a test of whether or not a gas 
cloud can provide pressure support against gravitational col- 
lapse. If a sound wave cannot travel across the cloud in the 
time it would take the cloud to gravit ationa l ly free fall to the 
centre, then the cloud will collapse. iKat j l)l992f) proposed 
that such a criterion should take the form: 

- > -f^= (3) 

where d is the sound speed of the gas particle in question 
and G is the gravitational constant. This tests if the tem- 
perature and pressure of one particle inside its smoothing 
sphere would be able to support t he whole sphere against 
gravitational collapse. IKatzl lll992f) only applies this crite- 
rion wh en the region is not affe cted by gravitational soften- 
ing, but lOkamoto et al.l i2005t) notes that the Jeans Crite- 
rion fo rmulated in this way is similar to the lBate fc Burkerd 
( 1997) resolution limit for artificial fragmentation. In H2.ll 
it is shown that this formulation of the Jeans criterion intro- 
duces a resolution dependence and hence we will eliminate 
it from our final formulation. 



2.2 Stochastic Star Formation 

Ideally, stars should form whenever spe cific crit e ria like 
th ose described above are m e t (see also iLi et alJ ll2005all 
or iMichel-Dansac fc Wozniakl ^0041) 1. However, since our 
simulations of galaxies have limited resolution we can only 
capture the global behaviour of star formation, and we use 
a pr obabi l istic approach. 

iKatzl il992T) bases the number of stars t hat fo rm o n the 
theoretical star formation rates of iLarsonl (ll96STt and ISilkl 
il987h . who both proposed that psfr ~ Pgas , where p is the 
volume density. These formulations make use of the fact that 
the dynamical time, tdvn ~ P~ 1|/ 2 - However, to form stars 
gas must also cool so IKatzl Jl992l) chose the star formation 
timescale, iform, to be the maximum of the dynamical time 
and the cooling time. If the gas was already co ol enough to 
form stars, i.e. T < T max , then t dyn was used. lKat3 (ll992Ti 
took T max = 30, 000K, which was appropriate given his cool- 
ing function. So one can write the star formation rate as 
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Figure 1. The mean star formation rate (SFR) as a function of 
the star particle mass in units of the original gas particle mass. 
The SFR shows little dependence on the initial star particle mass. 
These results are for our fiducial values of parameters and criteria 
on the isolated model Milky Way (IMMW) described in J2| 



dpi, 
~~dt 



-. Pgas 
^form 



(4) 



where we have introduced a constant efficiency factor c* that 
will enable us to adjust the star formation rate to match 
observations. We will choose to set tf orm = tdyn at all times 
and instead will introduce another star formation criterion 
that T < T max . 

Because of the dependence on density, it is possible to 
create a stochastic recipe for when and where stars should 
form. If one takes the probability that a star will form as 



-c*At/t to „ 



(5) 



where m gas is the mass of the gas particle and m star is the 
mass of the potentially spawned star particle then on average 
one recovers Equation (QJ. 

In this formulation there is a greater probability of a 
star forming in denser areas. For each star formation eli- 
gible gas particle a random number, r, is drawn between 
zero a nd on e and if r < p a new star particle is created. 
IKatzl lll992l) assumed that the star particle mass was al- 
ways a constant fraction of the parent gas particle mass, i.e. 
m s t a r/m gas = e*. In that case, e* _1 replaces m gas /m s t a r in 
Equation JKJ. e* can be absorbed into the definition of c* , 
for the relevant case where At /t f orm is small , and hence 
does not appear i n Equation ( 5) of IKatzl il992l) . Therefore, 
to compare with IKatzl il992l) . c*'s here should be multi- 
plied by e*. Here we are interested in making the new star 
particle have a constant mass to remove the large variation 
in star particle mass that can occur as stars form at dif- 
ferent times during the simulation. Using Equation (JKJ and 
assuming a constant stellar particle mass we can recover the 
same global star formation rate, almost independent of that 
mass, as shown in Figure A smaller star particle mass 
gives us better resolution for the stellar component at the 
cost of greater computational expense and we choose a value 
of 0.2 x the initial gas particle mass as a compromise. 
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Cod es that use multiphase gas pa rticles for star forma- 
tion, e.g. ISpringel fc Hernauisd i2003t) 'explicit' star forma- 
tion, more finely time sample star formation since a spe- 
cific amount of gas is converted into stars during every gas 
particle timestep and new star particles are only spawned 
periodically. However, this approach has an inherent prob- 
lem. When feedback occurs, stars that are still tied to their 
parent gas particles can get hydrodynamically accelerated 
and possibly even be ejected from the galaxy, which is not 
physical. 

One problem that could arise when using a stochastic 
approach for star formation is that particularly dense regions 
might not experience the vig orous star formation that is ob- 
served in starburst galaxies jTelescolll983 IChapman et alJ 
2003) if the mass of all star particles is fixed. A possible 
solution to this problem is to make the mass of the star par- 
ticle that for ms dependent on proper t ies of the gas particle. 
For instance, lElmeereen fc Efremovl |l997) provide a for- 
mula using the density and the pressure of gas to determine 
what fraction of the gas will form stars. Initial tests using 
this formula suggest that our resolution is too coarse for 
this recipe to yield useful results. However, for our fiducial 
parameters, a star particle mass 0.2 times the initial gas par- 
ticle mass and a star formation interval, At, of 1 Myr, this 
problem would only arise if gas consumption times become 
less than 5 Myr, a situation that is rarely, if ever, observed. 



3 FEEDBACK 



iKatd lll992T) and KWH add mass and energy feedback from 
type II supernova. The energy is added gradually with an 
exponential decay rate of 20 Myr. It is added at the lo- 
cation of the parent gas particle and is smoothed using the 
SPH smoothing kernel. Since this thermal energy is typically 
added to very dense gas, it is quickly radiated away and has 
little effect on the evolution of the galaxy. We describe a 
stronger and perhaps more realisti c method for in cluding 
type II supernovae (SNII) feedback iGerritsenlll997f) . which 
we implement here. In addition we now also include feedback 
from type la supernovae and stellar winds from planetary 
nebulae and allow the metals produced in our stars to be 
distributed. Their implementation is described towards the 
end of this section. 



3.1 Type II Supernovae (SNII) 

SNII play a dominant role in regulating star formation in our 
simulations because of their ability to heat large volumes of 
the interstellar medium near the site of star formation and 



consequently prevent more gas from collapsing ( Silk 2003 
Because the blastwave shocks of SNII convert the kinetic 
energy of ejecta into thermal energy on scales smaller than 
our simulations resolve, feedback in our simulations is purely 
thermal. 

The number of supernovae produced by a star particle 
depends on the initial mass function of the stars that form. 
We use the three pie ce power law fit of the IMF defined in 
iMiller fc Scale! ( 1979), where a = -0.4 for stars with masses 
between 0.1 and 1 Mq; a = -1.5 for stars with masses be- 
tween 1 and 10 Mq; a = -2.3 for stars with ma sses greater 
than 10 M e . Our IMF starts at 0.1 Mq because ItTeid et alJ 



(2002) report that the stellar luminosity function appears 
to turn over at a luminosity corresponding to that mass, 
so that stars with a mass < 0.1 Mq do not make a large 
contribution to the total stellar mass. 

The time when a sta r explodes as a supe rnova depends 
on its lifetime. We use the|Rai1^rj_et_alJ 1L996J) parameterisa- 
tion of the Padova group's iAlongi et alll993l : iBressan et alJ 
Il993t iBertelli et all Il994l) stellar lifetime calculations for 
stars of varying metallicities. In this parameterisation, since 
more massive stars have shorter lifetimes than small mass 
stars, it is possible to determine the maximum and mini- 
mum stellar mass that will explode during a given timestep 
and, therefore, integrating over the initial mass function 
provides the to tal mass and n umbe r of stars that will ex- 
plode. Like the iRaiteri et al] (I1996T) recipe, we only allow 
stars between 8 and 40 Mq to explode as SNII; stars more 
massive than this are assumed to either collapse into black 
holes or explode as type lb supernovae. Regardless, few stars 
form with masses greater than 40 Mq , so the impact on the 
feedback is minimal. The use of stellar lifetimes is an im- 
provement on previous SN feedback recipes that bled the 
supernova energy out gradually as so me type of exponential 
function after a star particle formed JCen fc OstrikeJll99A 
KWH). Thus, ma ny implementa t ions of feedback choose 
to use them now (iLia et al.ll2002l: iKawata fc Gibson! 120031: 
lOkamoto et aljl2003tlScannapieco et aljl200Sh . 

We multiply the number of SNII that explode by the 
energy ejected into the ISM, Esn, a fixed fraction of the 
canonical 10 51 ergs produced by a supernovae, and distribute 
that energy to the surrounding gas particles. The energy is 
spread out, weighted by the mass of the gas particle that 
receiv es it, using the SPH smoothing kernel. Unlike in lKatd 
( 1992) and KWH, however, the feedback is centred on the 
current position of the star particle and not on the parent 
gas particle. Therefore, the feedback energy received by a 
neighbouring gas particle i for a total feedback energy of 
A_Esn is just 

miW(\vi -r 8 |,A s )AJSsN 



(6) 



where h s is the distance from the star particle to its 32nd 
closest neighbouring gas particle. We explore how star for- 
mation depends on _Esn in When Esn = 10 50 ergs, 
7.65 x 10 47 ergs of energy are deposited into the surround- 
ing gas for every one Mq, o f stars formed. Metal production 
follows lilaiteri et alJ (I1996T) . where the ejected oxygen and 
iron mass are estimated using a power la w based on the mass 
that Raite rTet alJ (I1996T) estimated from lWooslev fc Weaver] 
(1995). We integrate this power law between the mini- 
mum and maximum mass stars that explode during a given 
timestep to determine the mass of oxygen and iron ejected 
into the ISM. The metals and the mass returned to the gas 
particles by the type II supernova are distributed in a way 
similar to Equation (|SJ. 

Unfortunately, the distributed feedback energy will have 
little impact on our simulated galaxies owing to our finite 
resolution and inability to resolve the complex multiphase 
interstellar medium properly. Since the supernovae explode 
in regions of high average density, the gas can radiate away 
the energy in much less time than a typical timestep. To 
make the feedback more realistic we tried two different ap- 
proaches (see H3.1.1l and H3.1.2H to mimic blastwaves, which 
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should take hundreds of thousands of years to cool. In both 
approaches we disable the radiative cooling in a number of 
the nearest gas particles. 

Disabling the radiative cooling in the surrounding gas 
particles allows us to model two different aspects of the feed- 
back phenomenon. First, the gas particles immediately sur- 
rounding new stars have their cooling disabled, and with the 
added energy from the supernovae, will likely become hotter 
than the maximum temperature, Tmax, allowed for forming 
stars. In this way, feedback inhibits further star formation 
in dense regions much like supernovae generate turbulence 
in molecular clouds, which provides global stability against 
further collapse of the entire molecular cloud. Secondly, the 
increased temperature of the gas particle models the high 
pressure of a blastwave, which plays a key role in shaping 
the interstellar medium, allowing the surrounding gas to nat- 
urally flow outwards. 

In our initial scheme, we determine the number of gas 
particles that have their cooling disabled by multiplying the 
SNII mass by a mass loading factor, j3, and shutoff the cool- 
ing for a constant time, r cso . Our other scheme requires fewer 
parameters because it depends on the analytic tre atment of 
blastwaves described in lMcKee fc Ostrikeil dl977f> . 

Both schemes allow gas particles to receive energy from 
multiple supernovae explosions in a similar fashion. Every 
timestep, each gas particle has its cooling shutoff time, rcso 
calculated from the total supernova energy received. For 
subsequent supernova, rcso is recalculated and extended 
when necessary. 



3.1.1 Supernova Mass Factor (f3) Recipe 

We derive t he mass loadi n g fact or co ncept from multi-phase 
recipe s like lYeoes et, alJ l)l997h and iHultman fc Pharasvnl 

( 1999) . We calculate the exact number of gas particles where 
we turn off the radiative cooling in the following manner. For 
gas particles surrounding a star particle, radiative cooling is 
disabled for gas particles within a sphere of mass /SMsnii- 
For a gas particle i, if 

/?MSNII > fa(|r ' 3 " rs|)3 pavc (7) 

where Msnii is the mass of stars that go supernovae during 
a given timestep, (|rj — r s |) is the distance from the star to 
the gas particle in question, and /9 avc is given by 

N 

Pave =y^TOjiy(|r s -rj|, ft,). (8) 
j=l 

The maximum number of particles for which we can disable 
the cooling is ATsmooth, which is 32 in our simulations. 

In our initial recipe, cooling is disabled for a fixed 
amount of time , rcso- We starte d wit h tcsq= 30 Myr based 
on the work of iGerritsenl il997f> and lThacker fe Couchmanl 

(2000) who suggest that 8 Mq stars have a lifetime of 30 
Myr (although using the Padova tracks it would be closer to 
38 Myr). Thus, after 30 Myr, feedback produced in a star 
forming region should be finished. We explore the effects of 
varying rcso in H5.2.2I 



3.1.2 Blastwave Recipe 

Exploration of the j3 parameter moti vated us to intro duce 
an explicit blastwave solu tion based on lChevalieJ (Il974l) and 
iMcKee fc Ostrikeil lll977l) . This solution reduces the num- 
ber of tunable parameters by providing both the maximum 
radius to which the blastwave explosion will reach and the 
time that the blastwave will keep the g as hot. The maxi mum 
radius of a supernova blastwave in the lChevalieJ Jl974f) sim- 
ulations was 

T> 1rl l. 74 t-,0.32 - 0.16 5-0.20 / n \ 

Re = 10 E 51 n P 04 pc (9) 

where Psn = -E51 10 ergs, no is the ambient Hydrogen den- 
sity, P04 = 10 _4 Pofc _1 where Po is the ambient pressure and 
k is the Boltzmann constant. Both no and Po are calculated 
using the SPH kernel for the gas particles surrounding the 
star. We temporarily turn off the radiative cooling for all gas 
particles within Re - However, there is an artificial maximum 
of 32 particles that can have their cooling disabled. Figure 
1181 shows how many of these particles there are typically. 

The simulations also provide a timescale for the time 
that a gas particle does not radiatively cool. Naively, one 
might think that this timescale should be the length of the 
Sedov phase of a supernova explosion. During this phase, en- 
ergy is conserved in the supernova remnant because it is not 
able to effectively radiate. However, the Sedov phase onl y 
lasts for tens of thousands of years dPadmanabhanl 1200 if) . 
Our simulations cannot resolve this timescale and the feed- 
back would be ineffective e ven if our simulations pr oduce 
clusters of supernova. Also, IMcKee fc Ostrikeil jl977l) sug- 
gest that a hot, low density shell survives well after the Sedov 
phase. 

Following the Sedov phase, blastwaves enter the snow- 
plow phase. During the snowplow phase, momentum is con- 
served as the blastwave expands b ecause the gas has cooled 
enough to radiate more efficiently. IMcKee fc Ostrikeil (Il977h 
present the end of the snowplow phase when a supernova 
remnant first reaches its maximum extent: 

, -1 n 5.92 r-,0.31 0.27 5-0.64 /, n N 

t E = 10 E 51 n P 04 yr (10) 

The supernova remnant continues to cool radiative ly even 
after it stops expanding. IMcKee fc Ostrikeil lll977f) report 
that the time that the hot, low density shell will survive is 

, , n 6. 85 t-,0.32 0.345-0.70 / 1 -1 \ 

tmax = 10 P51 n P 04 yr (11) 

Either of these timescales may be appropriate for the length 
of time to disable cooling and we report on how using either 
t e or tmax affects star formation in H5.3.1I 

Even in the case where a supernova has exploded in a 
previous timestep, we currently use the previous equations 
as described. We are investigating the possible interactions 
of multiple supernovae remnants. As an initial assumption, 
we suppose that all of the supernovae exploding during a 
given timestep combine their energy to generate the blast- 
wave. 

3.1.3 Small SN Smoothing 

Since we only disable cooling for a fraction of the particles 
within the smoothing radius, it is only those particles that 
maintain the high temperature generated from the super- 
nova. Thus, all the energy that gets distributed beyond the 
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blast radius is quickly radiated away, which is still unphys- 
ical. To address this problem, we introduce another variant 
of the blastwave approach where we restrict the distribu- 
tion of energy from the supernova only to those particles 
within the blast radius using a kernel function limited to 
just the particles within the blast radius. Initial trials only 
distributed metals and mass inside the blast radius like the 
energy. However, as we have yet to implement diffusion of 
metals between gas particles, the supernova explosion repre- 
sents the only time when metals can be widely distributed. 
Distributing metals only inside the blast radius lead to spu- 
rious metal distributions. Thus, we reverted to distributing 
the metals and mass across the entire smoothing sphere. 

It sometimes occurs that no particle is within the blast 
radius. In this case, we deposit the energy, metals, and mass 
to the nearest gas particle. Ejecta are distributed in this 
manner for both SNII and SNIa, but for SNIa, the cooling 
is not disabled. If there are no supernova ejecta, the wind 
feedback is distributed across all 32 nearest neighbours with 
the standard smoothing radius. 

Such an approa ch might be a cause for concern since 
IBenz fc Thielemann! jl99ftl have found that depositing en- 
ergy into a single gas particle in a SPH simulation can 
lead to overcooling and a violation of energy conservation 
owing to the large energy gradient introduced. However, 
ISpringel fc Hernauistl (I2002T) show that this is not a prob- 
lem if one uses the asymmetric form of the thermal energy 
equation, as we do in GASOLINE. 

3.2 Type la Supernovae (SN la) 

SNIa are also significant sources of metals and are thought 
to occur in binary systems. The method we use to deter- 
m ine how many SNIa explode is again described in detail 
in lRaiteri et al] (ll996T> . The minimum mass of a binary sys- 
tem is 3 Mq and the maximum mass i s 16 Mq (two 8 Mm 
stars) . Using the binary fractions from Rai teri et al.l l) 19961) 
makes the number of SNIa 10-20 % of the total supernovae 
in our Isolated Model Milky Way si mulations, as observed 
bv Ivan den Bereh fc McClurel (|l994t) in spiral galaxies. 

Usually we distribute the SNIa energy using the 
smoothing kernel amongst the JV smoot h nearest neighbour- 
ing gas particles. However, for those star particles that also 
distribute SNII energy within the blast radius, the SNIa en- 
ergy is also only distributed within the blast radius. 

Radiative cooling is not disabled as a result of SN la be- 
cause SN la occur much after the stars initially form. During 
this time, stars would dynamically spread out of their initial 
associations and subsequent SNIa would not be a collective 
phenomenon, like SNII, and hence would not lead to large 
blastwaves. Since stars in our simulations consist of indi- 
visible particles, the SNIa stars cannot drift apart as they 
should, would act collectively and produce too large an effect 
if we turned off the cooling. Our early simulations indicated 
the inclusion of SNIa produced too large of a feedback effect. 

Like energy, mass and metals are smoothed across all 
-^smooth nearest neighbour particles. All SNIa are assumed 
to eject the same mass (1.4 Mq) and the same amount 
of iron (0.63 Mm ) and oxygen (0.13 Mq) based on the 
iThielemann et all lll986t) SNIa yield models. These quan- 
tities are added to the existing oxygen and iron in the gas 
particles by mass and then converted to a fractional mass 



of the gas particle so that when a new star forms, it will 
form with the same fractional metal content as its parent 
gas particle. 

3.3 Stellar wind feedback 

The feedback contribution of stellar winds is also signifi- 
cant. Stars with masses below ~8 Mq return substantial 
fractions of their mass to the ISM as they evolve and leave 
behind white dwarf remnants. We ba se our wind feedback 
on the work of IKennicutt et al] (I1994T) who find that the to- 
tal stellar return fraction is 0.25 to 0.50 of the initial mass 
depending on the IMF. Because the return rate is so high, 
this form of feedback can greatly prolong star formation in 
galaxies without gas inflow. 

For simplicity, we consider only stars between 1 and 8 
Mq and assume that lower mass stars remain unevolved. To 
determine the fraction of mass returned for a given stellar 
mass we use the initial- final mass relation of IWeidemannI 
( 1987§) and then fit his results to a continuous function. 

In practical terms, we implement this feedback mecha- 
nism by first taking each star particle and determining the 
range of stellar masses t hat die during the cu rrent timestep 
using the lifetimes from lRaiteri et al] dl996i) . Then we cal- 
culate a returned mass fr action for th i s mas s range using 
the function derived from IWeidemannI il987T) . We add the 
feedback to the gas particles in the same manner as the SNe 
feedbacks, except without injecting any energy. The metal- 
licity of the returned gas is simply the metallicity of the star 
particle. In the future, we plan to include metal production 
by intermediate mass stars. The total fraction of mass lost 
from a star particle over >10 Gyr is 40% and of this ~99% 
of the mass loss results from stellar winds. 



4 TESTS OF THE RECIPE 

To closely examine our star formation formalism, we created 
an isolated model Milky Way (here afte r IMMW) based on 
the dynamical model presented in iKlypin et al] (l2002f) at 
several different resolutions. We use these models to tune 
our star formation recipe to produce results con s istent with 
observations (e. g. iRoc ha-Pinto fc Maciel Il997t iKennicutil 
ll998tlWong fc Blitzll2002D . 

4.1 Isolated Galaxy 

The IMMW was created using the specificati ons of |Sgringe 
(2000) in that it resides in a slightly modified iNavarro et al 
(1997) (hereafter NFW) dark matter halo where the cen- 
tral dark matter has been concentrated by infalling baryons. 
Thus, the density distribution and potential is slightly dif- 
ferent from a pure NFW halo. We initially modelled the dark 
matter using a velocity ^200 = 150 km s _1 , a concentration, 
c = r-2oo/r s = 12, resulting in a mass M200 = 1.12 x 10 12 M Q 
and radius r2oo=214 kpc at an average overdensity — = 
200. However, to r eplicate as close ly as possible the poten- 
tial created in the ISpringell (1200011 model, we fit a slightly 
different NFW profile to the particle potential. This fixed 
potential has M200 = 7.314 x IO^Mq, r 2 oo=153 kpc, and 
concentration c=20.648. 

The baryons are distributed in a stellar and gaseous 
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disk along with a spherical bulge that contains only star 
particles with a total baryonic mass of 4.55xlO lo M . The 
bulge contains 4.93xlO 9 M0, ~ 10% the baryonic mass. The 
stellar disk follows an exponential profile with a scale length 
of 3.5 kpc and constitutes 90% of the total mass of the disk. 
We distributed the remaining 10% of the baryonic disk mass 
as collisional SPH gas particles in an exponential profile with 
a sca le length of 7 kpc, tw ice the scale length of the stellar 
disk feroeils fc Rheelll997D . We introduce no holes or gaps 
in the gas because molecular gas studies show that where 
neutral atomic hydrogen column densities decrease in the 
inner par ts of disks, molecula r gas densities increase to fill 
the void dWong fc Blitzll2002l) . 

We evolve the disk until the initial instabilities (caused 
by the initially smooth particle distribution developing a 
spiral pattern) die out before turning on star formation. The 
minimum Toomre Q value is 2, so this disk is a very stable 
disk and grows no bar. 

To check whether our star formation rate converges to 
a single value at various resolutions, we replicate our model 
three times. The lowest resolution galaxy starts with only 
9,000 star particles of 5xlO 6 M and 2,000 gas particles of 
2.5xlO 6 M , the medium resolution galaxy has 45,000 star 
particles of lx 10 6 M Q and 10,000 gas particles of 5x 10 s M , 
and the highest resolution galaxy has 225,000 star parti- 
cles of 2xl0 5 M Q and 50,000 gas particles of lxlO 5 M . We 
soften the gravity using spline softening and our gravita- 
tional softening length, 650 pc for the gas particles and 
325 pc for star particles, remains fixed for all three resolu- 
tions. The equivalent Plummer softening is about 0.7 times 
smaller. 

4.2 Goals of IMMW 

Many observations have been dedicated to studying star for- 
mation in the Milky Way and spiral galaxies similar to our 
IMMW. The observations allow us to c onstrain the val ues of 
our star formation recipe parameters. ISchmidtl (ll95STl pro- 
vides the foremost constraint for how many stars should 
form in a dense gas environment. Schmidt showed that the 
surface density of star formation, Esfr, follows a power law 
of the gas surface den sity, £ gas , called t he Schmidt law. The 
more recent work of iKennicuttl l|l998f) specifies the exact 
slope and norma lisation of this relationship. Equation 4 of 
IKennicuttl \ 19981) states that: 

E SFR = (2.5±O.7)xlO- 4 (-|^) 1 ' 4±o - 15 M yr- 1 kpc- 2 (12) 
IMqPC A 

The formulation of our star formation recipe should en- 
sure that our star formation approximately follows the slope 
of this relationship, while the star formation efficiency, c*, 
should adjust the normalisation. 

Another constraint is the observed nearly steady star 
formation rate in the Milky Way. The local stellar neigh- 
bourhood shows evidence t hat the star formation rate 
has been constant for Gyr jRocha- Pinto fc Maciell Il997b 
when averaged over long timescales. In our simulations, a 
steady star formation rate results from the maintenance 
of a constant expon ential gas surface density profile. The 
IWong fc Blitzl l)2002l ) observations of exponential gas sur- 
face density profiles, therefore, provide another constraint 
related to the steady star formation rate. 



Our experiments consist of varying the four star for- 
mation criteria (temperature, density, converging flow, and 
Jeans) and the four parameters (c*, 15s n, A an< A r cso) to de- 
termine how each criterion and each parameter affects star 
formation. Criteria are solid cutoffs that eliminate gas parti- 
cles from forming stars whereas parameters are proportional 
constants that affect the rate of star formation and feedback. 
In this spirit, we choose a fiducial set of criteria and param- 
eters, and then proceed to vary each parameter or criterion 
individually. Our fiducial criteria are T max = 30,000 K, n m i n 
= 0.1 cm -3 , flows must be converging, and no Jeans-like 
criterion. The fiducial parameters are /3=10,000, c* = 0.1, 
tcso = 3 x 10 7 yr, and _Esn = 10 50 ergs. These parameters 
do not necessarily represent a best fit, which is presented in 
the conclusions, but are simply a starting point that produce 
relatively normal results. 



4.3 Numerical Precision 

Our star formation recipe i s implemented in th e parallel 
tree SPH code GASOLINE JWadslev et alJl2004l) . GASO- 
LINE implements cooling similar to what is described in 
KWH. It assumes ionisation equilibrium, an ideal gas with 
primordial composition, and solves for the abundances of 
each ion species. The scheme uses the collisional ionisation 
rates reported in lAbel etalj <TT997f). the radiative recombina - 
tion rates from Black! il98lf) and lVernerfc Fg^landl (1996), 
bremsstrahlung, and line cooling from ICenl (1992). The en- 
ergy integration uses a semi-implicit stiff integrator indepen- 
dently for each particle with the compressive heating and 
density (i.e. terms dependent on other particles) assumed to 
be constant over the timestep. 

Gravity is calculated for each particle using tree ele- 
ments that span at most 8 = 0.7 of the size of the tree el- 
ement's distance from the particle. Every particle has its 
forces calculated on each large time-step, 1.53 xlO 7 yr. 
GASOLINE is multistepping so that every particle's time- 
step Ai grav = Tj^J^, where r\ = 0.175, e< is the gravitational 
softening length, and <n is the acceleration. For gas parti- 
cles, the time-step must also be less than Ai gaa = ?7cour an t - 1 , 
where jycourant = 0.4 and c, is the sound speed. We restricted 
the smallest SPH smoothing length to be 0.01 of the gravi- 
tational softening length. 

Stars are formed and feedback is calculated every 1 Myr 
in the simulations. The time between star formation events 
has no relationship to the major timesteps of the simulation 
when every particle has its forces calculated. However, star 
formation is tied to the minor timesteps when some subset of 
the particles have their forces calculated. As these timesteps 
may not be exactly 1 Myr, we choose the minor timestep 
that is closest to the time that we want to form stars and 
add feedback at that time. 



5 RESULTS 

We have conducted a variety of experiments with our 
IMMW, adjusting the parameters discussed above. For each 
of these experiments, we have data on the star formation 
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Figure 2. A demonstration of unacceptable (top) and acceptable 
(bottom) star formation conditions. The top panels are the results 
produced when using a small star formation efficiency, c* of 0.01 
in contrast to the bottom panels that use the best fit value of 
c*=0.05 in the blastwave model. The left panels plot the SFR vs. 
time and the right panels the SFR surface density vs. gas surface 
density relation. The two dashed lines in the left panel show the 
region that is used to calculate the average SF Rs. The solid lin e 
in the right panel is the observed Schmidt Law <KennicutJl98dl . 



rate (SFR) as a function of time and gas surface density 1 . 
In Figure [5] we present two parameter choices and show 
star formation histories and SFR surface density versus gas 
surface density relations as an example of star formation 
variation within our experiments. In this case, the blast- 
wave feedback model is used in medium resolution (45,000 
stars, 10,000 gas) IMMWs. The top panels are the results 
produced when using a small star formation efficiency, c* 
of 0.01 in contrast to the best fit for the blastwave model 
c* =0.05 displayed in bottom plots. 

The star formation history, plotted in the left panels, 
is a histogram of when stars formed. The vertical dashed 
lines indicate the period between 6x10 s and 1.6xl0 9 yr af- 
ter the beginning of the simulation, the time over which 
we calculate the average SFR, which we use in subsequent 
plots. The simulations start 8 Gyr after the stars in the ini- 
tial conditions formed so that there are no feedback effects 
from those stars. We choose a time period past the begin- 
ning of the simulation because it falls well after any initial 
transient starburst that may result from non-equilibrium ini- 
tial conditions. Such a starburst can happen because all the 
gas particles at the start of the simulation are unaffected by 
feedback and hence all of them may be eligible to form stars. 
As shown in the panels, both choices result in an approx- 
imately constant star formation history but the amount of 
present day star formation in the upper panel is lower than 
the observed value. 

The right panels of Figure show SFR surface density 
versus gas surface density and indica tes how well our star 
formation follows a Schmidt law (IKennicutlill998r) . To cre- 
ate these plots, we azimuthally sum the mass of stars that 
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Figure 3. A comparison between star formation with (solid) and 
without (dotted) feedback. The run with feedback uses Esn = 
3x 10 50 ergs with the blastwave model at high resolution (275,000 
particles). 



formed over the last 100 Myr of the simulation in 500 pc 
radial bins and plot them against the final surface density 
of the gas in each bin. The panels show that our star forma- 
tion formulation creates a Schmidt Law with the right slope. 
However, the upper panel does not have the correct ampli- 
tude suggesting that the choice of parameters was incorrect. 
It is generally the case that an average star formation rate of 
around 0.8 Mq yr -1 in the high resolution case reproduces 
the Schmidt Law in the IMMW experiments. The turn down 
in SFR density at surface densities less than 2 p c~ 2 
matches the observations of iMartin fc Kennicutd j200ll) as 
we discuss in Section 

The effect of feedback is previewed in Figure [3] from 
a high resolution (225,000 stars, 50,000 gas) IMMW. The 
shallower drop-off in star formation rate with feedback owes 
in part to the fact that the supernova energy is being effec- 
tively used to prevent star formation and in part due to the 
stellar wind feedback that returns gas to the ISM. Both star 
formation rates decline exponentially because of the stochas- 
tic recipe ( H2.21 . Figure 2] shows how the feedback immedi- 
ately increases the temperature of the gas surrounding stars 
that have recently formed. The higher temperature results in 
higher pressure that allows the hot particle to push around 
the cooler gas surrounding it, as shown in Figure lT9l The ex- 
pansion leads to only a modest reduction in density because 
density is a smoothed property that includes the surround- 
ing high density particles. Thus, it is not low gas density 
that suppresses star formation, but high gas temperature. 

5.1 Effects of Criteria 

5.1.1 T max : Maximum Temperature 

As we al ready s tated , we added an additional criterion to 
those in iKata (|1992i) and KWH: gas particles may not 
form stars unless their temperature is below T max , typically 
10,000's K. This may seem like a high temperature threshold 
for star formation given that star forming molecular clouds 
are observed to cool down to ~ 100 K. However, our cooling 
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Figure 4. The phase diagram for the 50,000 gas particles in the 
high resolution (275,000 paticles) IMMW run with E SN = 3x 10 50 
ergs using the blastwave model. The fiducial values for T max and 
W m in are drawn to indicate the gas particles that pass the star 
formation criteria. All of the particles with temperatures above 
15,000 K are there as a result of feedback because all gas parti- 
cles in the disk start with T=10,000 K. A couple of gas particles 
show a modest decrease in density below the density threshold 
indicating that the gas does expand slightly owing to its high 
temperature and pressure. 
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Figure 5. The mean SFR as a function of T ma x at low resolution 
(asterisks), medium resolution (diamonds), and high resolution 
(triangles). 



is limited to H and He atomic cooling, which can only cool 
gas down to ~ 10, 000 K, and we average over scales much 
larger than star forming clouds. A future improvement to 
the code w il l be t o incl ude molecular hy drogen cooling (e.g. 
lAbel et alJ (Il997l) and iKravtsovl J2003h 'l. which will allow 
the gas to cool below 10,000 K, but even then, unless the 
resolution were greatly improved, T max should remain above 
10,000 K. 

Figure|^|shows the effect that varying T max from 10,000- 
50,000 K has on the mean SFR. As long as T max > 15, 000 



K the mean SFR is almost independent of T max . We might 
expect so little effect since most gas particles that are near 
star forming temperatures have cooled all the way to 10,000 
K. For example, only 18% of the gas mass in the medium 
resolution simulation, and only 20% of the gas mass at high 
resolution is warmer than 12,000 K after 1 Gyr of evolution 
using our fiducial recipe. It is only those gas particles that 
have been heated as the result of supernova feedback that 
are excluded from forming stars by the T max criterion even 
though they may remain in a dense environment. We tried 
lowering the threshold temperature all the way down to the 
mean gas particle temperature, 10,000 K. This produced 
galaxies with much burstier star formation histories. Gas 
particles would pile up just above the temperature thresh- 
old, cool all at once, produce lots of star formation, and then 
heat the gas up so that it could not form stars until the next 
cooling episode started the cycle all over again. This exper- 
iment demonstrates that T maj should stay above 12,000 K, 
but its specific value does not critically affect star forma- 
tion. We choose to use T max = 15, 000 K. We note that the 
T max criterion is critical to our feedback prescription because 
it enables star formation to be immediately suppressed by 
supernova feedback. 



5.1.2 n m i n : Minimum Density 

As n m i n increases, fewer gas particles are eligible to form 
stars, and hence fewer stars form. Figure |S| shows that a 
hundred-fold increase in the minimum density causes an or- 
der of magnitude reduction in the star formation rate. As 
the gas density profile declines exponentially, most of the 
gas particles that are not eligible to form stars reside in the 
outer parts of the disk. Because the IMMW disk is gravita- 
tionally stable, instabilities do not drive star formation and 
thus density simply correlates with radius. The minimum 
density that we choose, 0.1 cm -3 , confines star formation 
the inner 20 kpc of our model galaxy, which corresponds to 
surface densities of approximately 2Mq pc -2 and above and 
matches the observed minimum densi ty threshold for star 
forma tion in nearby spiral galaxies iMartin fc Kennicutd 
1200 ll). 



5.1.3 Jeans Criterion 

The Jeans Criterion f ii2.ll is a test of whether or not a 
gas cloud can provide enough pressure support to prevent 
gravitational collapse and is commonly used to test whether 
gas can form stars. Our initial recipe implemented it as a 
comparison between the sound crossing time, hi/ci, where 
hi is the smoothing length and a is the sound speed, and 
the dynamical time, l/^/inGpi. 

Figure [7| shows this comparison for all the gas particles 
in the simulation at three different resolutions. What is im- 
mediately apparent is that the sound crossing time decreases 
significantly with increasing resolution, which happens be- 
cause a region at given density is sampled by more particles 
in the higher resolution simulations. More resolution, with- 
out a corresponding drop in temperature, results in a smaller 
smoothing length making the sound crossing time shorter. 
This decrease in the sound crossing time becomes signifi- 
cant enough in the highest resolution simulation to make 
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Figure 6. The mean SFR as a function of minimum density 
(Pmin) a * low resolution (asterisks), medium resolution (dia- 
monds), and high resolution (triangles). 



gas particles unable to pass the Jeans Criterion, at least us- 
ing our formulation of it and the atomic transition cooling 
function employed in these simulations. We also plot a line 
representing particles at 9500 K in the high resolution sim- 
ulation, as cold as they can get using atomic cooling, using 
typical smoothing lengths, which shows that typical parti- 
cles have little chance of passing the Jeans Criterion unless 
they cool a great deal adiabatically. Gas at the highest densi- 
ties is sampled at the highest resolution but because the gas 
does not realistically cool, the sound crossing time becomes 
so short that td yn is still greater than hi/ci, and the dense 
gas unrealistically fails the Jeans Criterion. This might not 
be true if we included more realistic molecular cooling. 

We plot the star formation histories obtained when we 
include the Jeans Criterion in Figure [S] The consequences 
of this problem are evident as few stars are able to form 
in the highest resolution simulation. Thus, we reject using 
the Jeans Criterion and use the simpler collapse conditions 
imposed by requiring that the gas have a minimum density 
and a maximum temperature. 

Ideally the simulations would also include molecular 
cooling so that particles would have more realistic tempera- 
tures, making the Jeans Criterion a more appropriate test. 
Figure |7| shows how the cooling curve used in these simula- 
tions can only cool particles down to 9500 K at high reso- 
lution. Any further decrease in temperature is the result of 
adiabatic cooling. Figure [7] also shows how the temperature 
maximum plays a similar role to the Jeans criterion. How- 
ever, the limit changes with resolution in this plot because 
the sound crossing time decreases with decreasing smoothing 
length while the sound speed of 15,000 K gas remains con- 
stant. For clarity, we only plot the high resolution maximum 
temperature limit. In the current generation of simulations 
of forming galaxies, the models only attempt to replicate the 
generic properties of star formation that s eem to be well rep - 
resented by a Schmidt Law as observed bv lKennicutd lll99Sl) . 
A temperature maximum of T max =15,000 K ensures that 
only reasonably cool gas particles form stars. The stochastic 
implementation of the star formation formula ensures the 
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Figure 7. The sound crossing time across a smoothing length 
versus the local dynamical time for 3 different simulation resolu- 
tions. An increase in resolution, which decreases the smoothing 
lengths and hence the sound crossing times takes gas particles 
from satisfying the Jeans Criterion in low resolution simulations 
(blue points) to failing the criterion in high resolution runs (green 
points) because the atomic cooling included in these simulations 
only allows gas to cool to 10,000 K. The line for particles at 
9500 K uses typical smoothing lengths of 360 pc at 0.1 cm -3 
(*dyn ~ 80 Myr) and 92 pc at 6 cm -3 (i dyn = 10 Myr). We also 
plot the maximum temperature criterion for the high resolution 
simulation. 
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Figure 8. The mean star formation rates at three different resolu- 
tions when the Jeans criterion is enabled. While the low (dashed) 
and medium (dotted) resolution simulations form a similar num- 
ber of stars, star formation is nearly eliminated in the high reso- 
lution (solid) simulation except for a couple of bursts. 



correct slope and with n m i n = 0.1 cnrt 3 we reproduce the 
correct low density threshold, making a Jeans criterion un- 
necessary. 
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5.1.4 Converging flow 

Most of the particles that satisfy the temperature and den- 
sity criteria also satisfy the criterion that the flow be con- 
verging. Therefore, there is not a large difference in the num- 
ber of stars formed in our IMMW galaxy with the criterion 
turned on or off. The difference was minor enough that we 
decided to run all of the simulations reported in this paper 
with the converging flow criterion enabled since it might 
still be of benefit in nonequilibrium situations that we will 
explore in future work. 



5.1.5 Summary of Criteria 

The criteria that we settled on are: 



15,000 K. 

-3 



• The gas particle must be colder than T m 

• The gas particle must be denser than n m i n = 0.1 cm" 

• The gas particle must be overdense enough to be part 
of a virialised structure. 

• The particle must be part of a converging flow, i.e. V • 
v < 0. 

In addition we found that the Jeans Ins tability loses 
meaning at high resolution as it is formulated in lKata Jl992T) 
and that star formation is most sensitive to the n m i n crite- 
rion. 



5.2 Effects of Parameters in the /3 Model 

The criteria established in the above section are appropriate 
for either of our feedback schemes. In this section we discuss 
how parameter choice affects star formation in the f3 recipe 
( H3.1.H . For all of these simulations, we used our fiducial 
choices and changed only the specified parameter to measure 
its effect on the star formation. 



5.2.1 (5: Supernova Feedback Mass Factor 

The f3 parameter disti ngui shes our method from thos e 
used bvlGerritsenl (Il997l) and lThacker fc Couchmanl fcOOll) . 
iGerritserJ dl997l) turns off the radiative cooling only for one 
gas particle while iThacker fc Couchmanl i200 ill turn off the 
radiative cooling for all 32 particles within the smoothing 
radius. The /3 model gives us more flexibility. For example, 
turning off the cooling for all 32 neighbouring gas particles 
could produce a feedback that is too strong and resolution 
dependent. In a low resolution simulation, the smoothing 
length is much longer than in a high resolution simulation, 
so much more gas is affected and the feedback is stronger. 
However, this effect is mitigated to a certain extent since 
the star particles are more massive, so that the feedback is 
the result of more supernovae explosions. Another possible 
problem is that disabling the cooling for all 32 particles does 
not take into account how much energy has been released 
in the supernova explosions. Even if only one supernova ex- 
plodes, it has the same effect as if hundreds of supernova 
exploded. 

As expected, Figure[§]shows that the more gas particles 
that have their radiative cooling disabled, i.e. larger values 
of (3, the fewer stars form, because the gas particles are not 
able to satisfy the maximum temperature criterion. Once 
f3 becomes large enough such that all the 32 neighbouring 
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Figure 9. The star formation rate as a function of mass factor 
j3 at low resolution (asterisks), at medium resolution (diamonds), 
and at high resolution (triangles). 



gas particles that received supernova feedback energy have 
their cooling disabled, increasing f3 further has little addi- 
tional effect. The reduced star formation at very low j3 (the 
leftmost lowest resolution point in Figure [5J occurs because 
our analysis averages star formation after the initial burst 
where much of the gas is converted into stars leaving little 
gas to form stars at later times. 

Figure [TUl shows the details of the /3 feedback model and 
can be compared with similar plots for the blastwave model 
shown in Figure 1181 The plots show the results using the 
fiducial recipe with /3 = 10,000. The critical panel is in the 
upper, right-hand corner where we plot how many gas parti- 
cles have their cooling disabled for each star particle during 
a given star formation event. One can see that the number of 
particles affected is not resolution dependent. Star particles 
produce an amount of SN energy proportional to their mass 
(as in the upper- left panel of Figure IT^ and, therefore, one 
might expect the radii within which we disable the cooling 
to be larger for lower resolution simulations. However, the 
inter-particle separation also changes with resolution and 
thus the blastwave regions are proportionally the same size 
as shown in the lower left hand panel. Here we normalise 
the radii to the high resolution simulation by scaling with 
the mass in supernova to the one third power. To recover 
the actual radii, in the medium resolution simulation just 
multiply by 1.6 and in the low resolution simulation by 2.7. 
But since the (3 model does not take into account how blast- 
waves expand differently in different density environments, 
the blastwave model provides a more physical representation 
of supernova explosions as described in H3.1.2I 



5.2.2 t cso •' Cooling Shutoff Time 

We choose our fiducial rcso=30 Myr, the time that we dis- 
able the radiative cooling, ba sed UDon lOerritsenl l|1997l ) and 
IThacker fc CouchmarJ i2000f) . Figure fTTI shows that there is 
not a significant increase in the SFR as rcso decreases. We 
prefer a short rcso as it more closely matches the expected 
lifetime of SN blastwaves. 
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Figure 10. Detailed characteristics of supernova feedback in the 
/3 model for three different resolution simulations with the pa- 
rameters set at their fiducial values. Each resolution has widely 
varying particle numbers, so every histogram has been normal- 
ized to a maximum of 100. The statistics were compiled over the 
first 200 Myr of the simulation. The upper-left panel plots the 
distribution of energies that each star particle releases from su- 
pernova explosions every time supernova feedback is calculated (1 
Myr). That quantity is normalised using star particle rest mass 
energies (E| N /m*c 2 ) to make comparisons between resolutions 
easier. We note that the quantity Eg N , the energy the entire star 
particle releases into the ISM, is different than the quantity EgN 
used elsewhere in this paper for the energy that individual super- 
novae release into the ISM. Effectively, the normalised quantity 
represents the efficiency with which the matter in the star par- 
ticles are converted into energy. The upper-right panel plots the 
number of particles with their cooling disabled. The relationship 
between gas particle mass and resolution means that the num- 
ber of particles with cooling disabled roughly traces the gas mass 
with cooling disabled, which is similar for every resolution. The 
lower-left panel plots the distribution of radii within which the 
cooling is disabled renormalised to the high resolution simulation 
by scaling by the mass in supernova to the one third power. 
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Figure 11. The mean SFR as a function of the time that radia- 
tive cooling is disabled at low resolution (asterisks), at medium 
resolution (diamonds), and at high resolution (triangles). 
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5.2.3 c* : Star Formation Efficiency 

c* controls the efficiency with which stars form. The c* pa- 
rameter can be thou ght of as either modifying th e timescale 
for star formation dSpringel fc Hernauisd 12003) or as the 
fraction of gas that becomes stars. The distribution of star 
formation timescales, tt OTm , for gas particles that pass the 
star formation criteria is close to a log normal distribution 
with a peak around 20 Myr and a tail out to 80 Myr. The 
results are plotted for our fiducial values at medium reso- 
lution and the times range from 20 to 80 Myr. Remember 
that tform = tdyn = 1/ V in Gp so the times just trace the gas 
density. Our fiducial value of c*=0.1 either corresponds to 
star formation timescales of ~ 300 Myr or implies that 
of a gas particle gets converted to stars each star formation 
timescale, ~ 30Myr. 

Figure H2l shows that c* has the strongest effect on star 
formation of any of the parameters that we have investigated 
so far. The val ue for c* is thus tightly constrained using the 
Schmidt Law <lKennicutdll998() defined in Equation 1121 and 
the best fit for the /3 feedback model is with c* ~ 0.1. 



Figure 12. The mean SFR as a function of c* at low resolu- 
tion (asterisks), at medium resolution (diamonds), and at high 
resolution (triangles). 



5.2.4 Esn-' Supernova Energy 

During the late stages of its life, the core of a massive star 
possesses 10 53 ergs of gravitational potential energy. How 
much of that energy is converted into thermal energy in the 
surrounding interstellar medium during a Type II SN explo- 
sion remains an open question. Most of the energy is con- 
verted into neutrinos and a small fraction of the neutrino 
flux interacts with enough matter to pro vide the canoni- 
cal SN kinetic energy value of 10 ergs (IColgate fc White! 
1966). Integrating SN light curves reveals that only 10 49 erg s 
are radiated away in the initial explosion jFilippenkdll997fi . 
but it is not clear that the rest of the energy is transferred 
to the ISM. In the current recipe, we have chosen to leave 
the energy transfer efficien cy as a parameter, Esn, and use 
the lThornton et alJ <ll99Sf) estimate that 10% of the super- 
nova's kinetic energy is converted into thermal energy as 
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Figure 13. The mean SFR as a function of -Esn using the j3 
feedback model at low resolution (asterisks), at medium resolu- 
tion (diamonds), and at high resolution (triangles). 
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Figure 14. The mean SFR as a function of c* at three differ- 
ent resolutions using the blastwave model and smoothing over 32 
particles for simulations at low resolution (asterisks), at medium 
resolution (diamonds), and at high resolution (triangles). 



our fiducial value. The rest of the energy is radiated away 
throughout the lifetime of the supernova blastwave. 

We distribute the supernova energy across all 32 nearest 
neighbour particles using the smoothing kernel. Since we are 
only disabling the cooling in a limited set of those particles, 
some of the energy will be quickly radiated away and have 
no effect on the simulation. However, the majority of the 
energy goes to particles that do have their radiative cooling 
disabled and thus the feedback does effect star formation. 

Figure IT3l shows that more stars form when less super- 
nova energy is returned to the ISM, as one would expect. 
All of the Esn values came close to fitting the Schmidt Law, 
with 2 x 10 50 < E SN < 6 x 10 50 doing the best. E SN does not 
have its largest impact on star formation in our IMMW. Its 
value is better constrained using other observables (see H6.1H 
or smaller galaxies. The reason that Esn effects SFRs for the 
IMMW is that gas particles get heated more when Esn is 
increased and thus are less likely to fall below the T max . The 
gas may also expand more owing to the additional pressure 
and be less likely to satisfy the n m in criterion. We note that 
the star formation rate begins to converge to one mean SFR 
as -Esn approaches 10 51 ergs. 

5.3 Effects of Parameters on the Analytic 
Blastwave Model 

The results of the previous section proved that turning off 
the radiative cooling of gas particles is an effective means 
of implementing supernova feedback in our IMMW. The re- 
sults also show that shorter tcso and moderate j3 values 
are in good agreement with the analyt ic blastwave solutions 
presented in lMcKee fc Ostrikerl jl977l) . Using these analytic 
expressions (detailed in H3.1.2H leaves only two free param- 
eters, c* and .Esn, governing both star formation and feed- 
back and is well motivated physically. 

The results of two different versions of the blastwave 
model are presented below. First, supernova feedbacks, i.e. 
thermal energy, metals, and mass, are distributed across 
the entire smoothing radius as described in H5.3.1I However, 



since energy that is added to particles that can cool is lost 
almost immediately we propose a second variant where we 
concentrate the feedback to only those particles that have 
their cooling disabled as we describe in H5.3.2I 

5.3.1 Smoothing over 32 Particles 

The star formation efficiency, c* , has a large effect on star 
formation when we use the analytic blastwave model much 
like in the (3 model as we show in Figure im Most values of e* 
do not re produce the obse rved Schmidt Law; only c*=0.05 
gives the iKennicuttl dl99ctl normalisation at all three reso- 
lutions. Therefore, every simulation described hereafter uses 
c* = 0.05 as its fiducial value. 

As discussed in H3.1.2I there are two possible shutoff 
times in the blastwave method, tE, the time that the snow- 
plow phase ends, and i max the time that it takes the gas to 
cool down to the ambient ISM temperature. We find that in 
the IMMW simulations there is little difference in the SFR 
between the two different choices with each having a SFR 
of 0.8 Mq yr _1 at high resolution with c* = 0.05. 

The amount of supernova energy transferred to the ISM 
has little effect on star formation as shown in Figure All 
of the values of Esn using c*= 0.05 fit the observed Schmidt 
Law well. Supernova energy provides local feedback and this 
maintains a steady SFR, but the amount of energy does not 
have a huge impact on the amount of star formation. 

5.3.2 SN Ejecta Smoothed only over the Blast Radius 

In our effort to inject the energy from supernovae ejecta into 
the ISM both efficiently and more physically, we concentrate 
the energy into just those particles inside of the blast radius, 
i.e. those that will have their cooling disabled when SNII 
explode. We also concentrate the energy and metal deposi- 
tion for SNIa as a matter of convenience, but do not disable 
radiative cooling for gas particles that are only within the 
SNIa blast radius. When there are no gas particles inside the 
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Figure 15. The mean SFR as a function of supernova energy for 
the blastwave feedback model with smoothing over 32 particles 
(thick, blue) and when the feedback is concentrated (thin, red) at 
low resolution (asterisks), at medium resolution (diamonds), and 
at high resolution (triangles). 



blast radius, as often occurs when for stellar wind feedback, 
we deposit all of the energy into the nearest gas particle. 

Concentrating the ejecta inside the blast radius has an 
impact on star formation. The average star formation rate 
drops to 0.7 M Q yr" 1 at c*=0.05 and £? S n = 10 50 ergs for 
the high resolution simulations from the 0.8 Mq yr _1 shown 
previously. While the star formation decreases, the shape of 
the star formation history does not change. 

Figure 1151 shows how the mean SFR decreases more 
when the energy is concentrated within the blast radius as 
opposed to spreading it over all 32 neighbouring particles. 
However, the variation in mean SFR does not lead to any 
deviation from the Schmidt Law, so it is not possible to con- 
strain the value of _Bsn from star formation in the IMMW 
alone. 



5.4 Effects of Resolution 

One trend that is apparent from all the plots is that more 
stars form at higher resolution. With the fiducial recipe, 
there is a factor of two difference in SFRs between the low- 
est resolution simulation with 11,000 particles and the high- 
est resolution simulation with 275,000 particles. There is a 
smaller difference between the medium and high resolution 
simulations than there is between the low and medium res- 
olutions so perhaps the results are converging. But a factor 
of two remains a significant impediment to comparing sim- 
ulations with observations. 

The difference in star formation is caused in part by the 
higher densities that gas particles reach in higher resolution 
simulations. Denser gas forms more stars because of Equa- 
tion 21 The bottom panel of Figure [TBI shows that the most 
critical density range is between 0.5 and 1.5 cm -3 . Most 
of the star forming eligible gas mass falls in this range and 
there is a significant difference in the fraction of gas with 
these densities. In the high resolution simulation two times 
more gas has these densities than in the low resolution sim- 
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Figure 16. The top panel shows how the SFR formation rate 
as a function of time increases with increasing resolution for the 
fiducial parameter settings of the /3 model at three different reso- 
lutions. The bottom panel is a cumulative plot of gas mass frac- 
tion at various densities and shows why the SFR increases, as a 
greater fraction of gas is at star forming densities in the higher 
resolution simulations. 
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Figure 17. The mean SFR as a function of the timestep used 
for star formation, Atsp, for the fiducial blastwave model at low 
resolution (asterisks), at medium resolution (diamonds), and at 
high resolution (triangles). 



ulation. There are smaller fractions of gas at high density 
and Figure \Wl shows that the different resolutions vary ran- 
domly in their high density gas mass content. The top panel 
of Figure Hoi shows that these variations have little bearing 
on how many stars form. 

Time resolution can also affect the SFR; we calculate 
star formation every A£sf. Our value of AtsF = 1 Myr 
comes from considering the fact that O star lifetimes can be 
as short as 1 Myr. Figure H7I however, shows that the actual 
value of AtsF is not crucial until either it starts to approach 
the time over which we disable the cooling, calculated using 
the analytic blastwave solution, which is typically around 10 
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Figure 18. The same as Figure UHl for the blastwave feedback 
model where all the feedback is concentrated within the blast 
radius. There is one additional panel in the lower right where we 
plot the distribution of the cooling shutoff times. 




Myr, or the length of a major system timestep, 15.6 Myr in 
these simulations. As star formation and feedback are not 
computationally intensive, calculating it every Myr is prac- 
tical and will not adversely affect star formation. 



Figure 19. An example of the effect that supernova feedback has 
on the structure of our IMMW. The figure measures 24 kpc on a 
side. Notice the hot, yellow particles that have pushed the colder, 
more purple particles more in the low density outer regions of the 
galaxy and less in the dense, inner regions. 



5.4-1 Feedback Behaviour 

The feedback method also contributes to the resolution sen- 
sitivity. As in the (3 method, more particles have their cool- 
ing disabled when the blastwave recipe is used at low reso- 
lution than at higher resolution as shown in Figure HS1 The 
radial extent of the blastwave depends on the energy of the 
explosion, so the bigger star particles in the low resolution 
simulation generate a larger effect than the larger number of 
star particles in the higher resolution simulations. However, 
also like in the f3 method, the differences disappear if one 
rescales by the mass in supernova to the one third power. 
The length of time that we disable cooling does depend on 
the resolution, i.e. the cooling is disabled for longer peri- 
ods at lower resolution. This is because supernova feedback 
events are larger but less frequent at lower resolution. At any 
one time there is about the same amount of gas mass with its 
cooling disabled at all three resolutions. It should be noted 
that Figures l9l II II [T3l and 1 151 show that star formation con- 
verges as the amount of feedback increases. Unfortunately, 
these values for the feedback are unphysical and produce 
fewer stars than observed. 

6 DISCUSSION 
6.1 Constraints 

The iKennicuttl (|l99ST l and iMartin fc KennicuttJ l|200lll ob- 
servations provide three constraints on how the star for- 
mation rate relates to the gas surface density. The first 
is that the logarithmic slope is 1.4±0.2 and constrains us 
to use the stochastic formulation of star formation where 



SFR~ p/tdyn- The second is the normalisation of that rela- 
tion. Using the j3 recipe, there are many different parameter 
combinations that can produce the proper normalisation. 
Using the blastwave method, only c* gre atly influen c es the 
star formation rate and c*=0.05 fits the iKennicutd dl99g|) 
observations t he best. Finally, there is the low density cut- 
off observed bv lMartin fc Kennicutd ( 1200 ll) . which we repro- 
duce with our n m i n criter ion of 0.1cm~ 3 . 

The observations of iRocha-Pinto fc Maciell il997T) in- 
dicate that stars have formed at roughly a constant rate 
for some time in the Milky Way based on the metallicity 
of G dwarfs. In the isolated model Milky Way (IMMW) 
star formation stays relatively cons tant with a slow decay 
when feedback is enabled. Figure 3 of lRocha-Pinto fc Maciell 
( 1997) provides no conclusive evidence of whether or not the 
star formation rate declines significantly. The major differ- 
ence between our simulations and the Milky Way is that the 
IMMW is isolated. It has no external source of gas like the 
real Milky Way, which may be in the form of gas rich merg- 
ers or cold gas streaming in along filaments, and it is not 
gravitationally disturbed by the presence of orbiting satel- 
lites. 

Figure 1191 shows that our feedback method produces 
hot gas particles that push the cold gas into dense filaments 
where the gas particles are cool and dense enough to form 
stars. The holes that the SN blow open are smaller in the 
dense centre of our IMMW and grow progressively larger 
as the SN explode in the less dense regions in the outskirts 
of the galaxy. Images of the Large Mag ellanic Cloud in H I 
show similar features (se e Figu re 4 of iJones et alJ (Il999ll 
or Figure 3 of iKim et al.l <|2003T) for a similar i mage of the 
Circinus Galaxy). Inside of the Milky Way disk. iHartmannl 
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Figure 20. The radial velocity dispersion vs. radius of the gas 
particles in our concentrated blastwave model for four choices of 
the supernova feedback energy, Esn at the three resolutions. We 
calculate the dispersion within 10 evenly spaced, radial bins be- 
tween and 10 kpc. The dispersion is most accurately determined 
for the highest resolution simulation (solid line). 



( 2002) observes that star formation happens along filamen- 
tary structures. Because of the jostling provided by the adi- 
abatic expansion of supernova heated gas, the site of star 
formation is constantly shifting from one high density re- 
gion to the next high density region. Our simulations are 
not run with enough resolution to follow the fine turbulence 
and Shockwaves that lead to star formation, but do show 
structures that are reminiscent of what is observed. 

We cannot use star formation properties to constrain 
the one free parameter in the blastwave feedback method, 
the amount of energy per supernova, Esn, because in simula- 
tions of the IMMW the SFR is not very sensitive to changes 
in -Esn. To constrain this parameter we use the radial ve- 
locity dispersion of the gas in the disk plane and the mass 
fraction of gas in the hot phase. In Figurc l2*Ul we plot the ra- 
dial velocity dispersion of the gas versus radius for the three 
resolutions using four different values for Esn- The disper- 
sions are almost independent of radius as observed in spiral 
galaxies. They are higher for the lowest resolution simula- 
tion but similar for the two highest resolutions. As expected 
the dispersion increases with the energy of the feedback. 
Typical vertical velocity dispersions of gas (under the as- 
sumption that random gas motions are isotropic) measured 
in quiescent spiral galaxie s are a little more th an 10 km s _1 
JPickev fc LockmanllT990l: iDickev et alJll990D . so we prefer 
Esn ~ 10 50 . 

In Figure 1211 we plot the cumulative mass fraction of 
gas above a given temperature. We plot this at all three res- 
olutions for the same four choices of Esn as in Figure 1201 
The temperature mass fractions are almost independent of 
resolution except for the temperature of the hottest gas in 
the two cases with the smallest supernova feedback. In these 
two plots the higher the resolution the hotter the gas can be- 
come. More importantly, the fraction of gas in the cold star 
forming phase is very similar across the three resolutions. 
As expected, when the supernova feedback energy becomes 
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Figure 21. The cumulative mass fraction of gas above a given 
temperature for the blastwave method using four different super- 
nova feedback energies, Esn at each of the three resolutions. 



larger, the gas can become hotter and a smaller fraction of 
gas remains in the cold phase. Early soft x-ra y observations 
ijOorenstein et, alJll974l : iBurstein et al.l ll977l) and observa- 
tions of th e O VI absorption line Jjenkins fc Melovlir974l : 
lYorklll974l) placed limits on the mass fraction of hot gas in 
the Milky Way with a temperature greater t han 5 x 10 5 K. 
Hot gas mass makes up ~ 0.01 of the ISM llvan der Hulsd 
1996). Since the dispersion preferred value of Esn = 10 50 
agrees with these observations, we fix the supernova energy 
at this value. 



6.2 Live Halo 

All the simulations reported so far have used a rigid analytic 
model to represent the dark matter halo. Live dark matter 
halos could potentially introduce noise or allow secular in- 
stabilities to develop that could change our SFRs. To investi- 
gat e this pos s ibility we resimulate the IMMW created using 
the lSoringeJ ll200ff) initial conditions, but for these simula- 
tions we do not remove the dark matter. We use our fiducial 
concentrated blastwave feedback method and all the fiducial 
parameter values for the live halo simulations. We evolve all 
three resolutions: for the high resolution simulation we use 
one million dark matter particles to represent the dark halo 
and use the same relative number of particles for the other 
two resolutions. As we show in Figure 1221 the addition of 
a live halo does not significantly change the mean SFRs at 
any resolution. 



6.3 Summary of Chosen Parameters 

Based on the constraints discussed in H6.ll we determined 
the parameter values in our star formation/supernova feed- 
back r ecipe. We eliminated the Jeans Criterion used in iKatzl 
( 1992) because it was sensitive to variations in resolution. 
Ou r parameter study using t he f3 model converged towards 
the lMcKee fc Ostrikerl (|l977f) analytic blastwave solution, so 
we chose to use the blastwave feedback method and concen- 
trated the feedback to those particles that have their cooling 
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Figure 22. The mean SFR for the final blastwave model both 
with and without a live dark matter halo. 



disabled. Employing this method means that there are only 
two free parameters, c* and Esn- We constrained the star 
formation e fficiency, c* to 0.05, which is close to the effi- 
ciency that lLada fc Ladal (I2003T) find in molecular clouds 
and star clusters. The mass of gas particles in the IMMW 
are close to the mass of molecular clouds, so this value for 
c* seems reasonable. The amount of supernova energy trans- 
ferred to the ISM, Esn, is difficult to constrain using SFRs 
in a massive system like the IMMW since the galactic po- 
tential is so much deeper than the amount of energy that 
SN can inject into the ISM. The value of £ S n = 10 50 ergs 
comes from comparing gas velocity dispersions and temper- 
ature distributions with observed values. 

The blastwave model also requires that three criterion 
are met before stars can form: T max , n m i n and that the gas 
is in an overdense virialised region. Gas particles may only 
form stars when their temperatures are below 15,000 K and 
is critical for making our SN feedback mechanism effective. 
Gas particles must also have a density above n m i n =0.1 cm -3 . 
The density criterion limits star formation to the dense 
regions of galaxies, which corresponds well to the d ensity 
limits observed bv iMartin fc Kennicuttl <|200ll l. The iKatzl 
(1992) argument that gas particles be part of a converging 
flow remains in our recipe because it does not have a major 
effect on the number of stars that form and might be useful 
in nonequilibrium situations. 



(1997) resolution criteria. When simulations are run with 
an isothermal equation of state, the gas collapses because 
of instabilities into dense clumps where stars must form. 
ILi et al] ll200rf) find that the coll apse happens with the exact 
surface density d ependence that lKennicuttl (I1998T) observed. 
IKravtsovl l)2003|) notices a steepening of the star formation 
rate with density in his simulations. The steepening led him 
to use a constant star formation timesca le rather t han a 
dynamical time that depends on density. IKravtsovl (|2003l) 
argues that this type of behaviour is the result of a tur- 
bulent buildup of high density gas in such a way that the 
higher the density of t he gas, th e more high density, star 
forming gas is present. IKravtsovl (120031) only uses thermal 
energy from supernovae as feedback, and hence it has little 
effect in his simulations. His need for altering the amount 
of star formation arises simply from how the gas collapses, 
and his adaptive mesh code may be more effective at re- 
solving turbulence than our SPH code. However, the idea 
of changing the relationship between star formation and gas 
density could be useful for creating a recipe with realistic 
star formation at very high resolutions. 

ITasker fc Brvanl l|200fift run a set of simulations similar 
to lLiet£dT(l2005bl) except that they include supernova feed- 
back and the gas is simulated on a grid. Even though their 
feedback implementation simply pours the feedback energy 
into the ISM, it proves effective at limiting star formation. 
The effectiveness of the feedback may be the result of high 
resolution, effective shock capture, and a more realistic cool 
curve using the Enzo AMR grid pr ogram. Even at the high 
resolutions, ITasker fc BrvanT(|2006l) find it necessary to use 
a recipe like Equation [IJthat fixes the star formation depen- 
dence on density to a Schmidt law. 

Recent obs ervations also prov i de cl ues to a new 
physical recipe. iBlitz fc Rosolowskvl ((2004) base a recipe 
on molecular gas observatio ns. The recipe is much like 
lElmegreen fc Efremovl dl997h in that it identifies gas pres- 
sure as the property most responsible for predicting SFRs. 
Observations also show interesting star formation behaviour 
in merging galaxies, the place where the most vi g orous star 
formation happens in the local universe. iBarnesI l|2004l) has 
pointed out that current star formation recipes have a dif- 
ficult time creating the starbursts observed in the shocked 
regions of merging galaxies, since star formation remains fo- 
cused in the central regions of galaxies where the densities 
are high enough to facilitate star formation. It is yet to be 
seen whether or not our star formation formulation suffers 
from this problem. 



6.4 Comparison to Other Work 

Several authors have recently proposed alternati ves to the 
star fo rmat ion recipe describ ed in this paper. Both lKravtsovl 
(2003) and lLi et all l|2005al) have run simulations in which 
they do not impose a Schmidt Law for the star formation. In- 
stead, they obtain this result naturally. While their studies 
shed considerable light on how stars form in the environ- 
ment of a galaxy, our goal is to simply produce reasonable 
star formation rates in simulati ons and thus relies on the 
observati ons of | K ermigy£ttj l)l998l) . 

The ILi et all l)2005af) recipe relies on high resolution 
simulations that at minimum satisfy the iBate fc Burkerd 



7 CONCLUSIONS 

Effective comparison of simulations with observations re- 
quires that simulations include gas and stars. Computa- 
tional limitations prevent simulations from representing ev- 
ery atom or even every star in the uni verse. However, there 
are global properties of star formation (lKennicuttll998l) that 
can be reproduced in simulations . The simulat ions presented 
here use a scheme presented in IKatzl IJ1992T) as a starting 
point. The maximum temperature, minimum density, con- 
verging flow, a nd sto c hastic selection of star forming gas par- 
ticles from the IKatzl dl992f) model were retained in our star 
formation recipe. However, the Jeans Criterion presented 
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in that work showed a strong resolution dependence in our 
simulations and was eliminated from our final star formation 
method. 

Simulations of an isolated model Milky Way (IMMW) 
that only included this simple star formation recipe and 
pure thermal energy supernova feedback had SFRs that were 
higher than those observed. Therefore, we implemented an 
effective version of supern ova feedback that uses th e blast- 
wave solution presented in lMcKee fc Ostrikerl l|l977T) . Simu- 
lations cannot resolve the blastwave shocks from supernovae. 
Since star formation occurs in dense gas regions, when the 
supernova energy is simply distributed amongst nearby gas 
particles as thermal energy, it is quickly radiated away and 
produces no effective feedback. Thus, we chose to disable 
the radiative cooling of gas particles in the proximity of re- 
cently exploded supernov ae so that the supernovae wil l sup- 
press star formation. The Th acker fe Couchm an (2001) sim- 
ulations showed the promise of this method, and we added 
the flexibility of how many gas particles have their cool- 
ing disabled to elimin a te the resolution dependence of the 
iThacker fc Couchmar] l)200lfl method. 

After an initial attempt using a recipe containing many 
parameters (our j3 model), we se ttled on the treatment o f 
supernova blastwaves presented in lMcKee fc Ostrikerl lll977l) 
to determine how many particles would have their cooling 
disabled. This proved an effective means for regulating star 
formation, and we were able to reproduce the Schmidt Law 
observed bv lKennicuttl l|l99Sh with c*=0.05. Other than the 
feedback effect necessary to make the IMMW's star forma- 
tion constant, the supernova feedback has little impact on 
star formation and we expect it to have little effect in other 
massive systems. 

Our final star formation recipe is made up of these pa- 
rameter values: 

• Tmax (maximum temperature) = 15,000 K 

• must be in a virialised region 

• iimin (minimum density) = 0.1 cm -3 

• must be in a converging flow, i.e. V ■ v < 

• -Esn (energy transfered from SN to ISM) = 10 50 ergs 

• c* (star formation efficiency) =0.05 

• R E (SN blast radius) = lO 1 ' 74 ^ 32 ^ 016 P^°- W pc 

• tmax (blast radius cooling time) = 

1n 6.85p0.32 0.34 p-0.70 

Given the success of our star formation method in the 
isolated Milky Way galaxy we hope it will continue in other 
situations. In future work we will investigate its performance 
when the galaxies have lower mass and see whether or not 
our feedback prescription can blow galactic winds like those 
observed. We are testing it on much more poorly resolved 
galaxies, like those that commonly arise in cosmological sim- 
ulations. Finally, we are applying this method to cosmologi- 
cal simulations to see if we can reproduce the observed shal- 
low faint end slope of the galaxy luminosity function. 
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